Direct MD simulation of liquid-solid phase equilibria for two-component plasmas 
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We determine the liquid-solid phase diagram for carbon-oxygen and oxygen-selenium plasma 
mixtures using two-phase MD simulations. We identified liquid, solid, and interface regions using 
a bond angle metric. To study finite size effects, we perform 27648 and 55296 ion simulations. To 
help monitor non-equilibrium effects, we calculate diffusion constants Di. For the carbon-oxygen 
system we find that Do for oxygen ions in the solid is much smaller than Dc for carbon ions and 
that both diffusion constants are 80 or more times smaller than diffusion constants in the liquid 
phase. There is excellent agreement between our carbon-oxygen phase diagram and that predicted 
by Medin and Gumming. This suggests that errors from finite size and non-equilibrium effects are 
small and that the carbon-oxygen phase diagram is now accurately known. The oxygen-selenium 
system is a simple two-component model for more complex rapid proton capture nucleosynthesis 
ash compositions for an accreting neutron star. Diffusion of oxygen, in a predominately selenium 
crystal, is remarkably fast, comparable to diffusion in the liquid phase. We find a somewhat lower 
melting temperature for the oxygen-selenium system than that predicted by Medin and Gumming. 
This is probably because of electron screening effects. 

PACS numbers: 97.20.Rp , 64.70.D- , 64.70.dg 



I. INTRODUCTION 

Observations of cooling White Dwarf (WD) stars pro- 
vide important information on the ages and evolution of 
stellar systems W-A . The interior of a WD is a Coulomb 
plasma of ions and a degenerate electron gas. As the 
star cools this plasma crystallizes. This can delay WD 
cooling, see for example ref. 0. Winget et al. recently 
observed effects from the latent heat of crystallization on 
the luminosity function of WDs in the globular cluster 
NGC 6397 jQ. Winget et al.'s observations may con- 
strain the melting temperature of the carbon and oxygen 
mixtures expected in these WD cores. In addition, as- 
troseismology provides an alternative way to study crys- 
tallization in WD, see for example fP. 

Furthermore, material accreting onto a neutron star 
(NS) will freeze to form new NS crust. A variety of 
nuclear reactions can take place, including rapid proton 
capture nucleosynthesis (rp process [9]) followed by 
electron captures, as the material is advected to higher 
densities. Horowitz et al. studied the crystallization of a 
complex rp process ash consisting of 17 chemical elements 
from oxygen to selenium [10 . They found chemical sep- 
aration upon freezing, with low Z elements preferentially 
remaining in the liquid NS ocean while high Z elements 
crystallize to form new NS crust. This change in com- 
position of the ocean may be important for superbursts. 
These are thought to be energetic thermonuclear explo- 
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sions of carbon pTMl3| . 

Melting temperatures and other properties of the 
liquid-solid phase diagram for multi-component plasmas 
have been determined from computer simulations. Segre- 
tain et al. calculated the carbon-oxygen phase diagram 
assuming a local density model for the free energy of 
the solid [H]. While, Ogata et al. [T1],[TS] and DeWitt 
et al. [17], [15] calculated the phase diagram based on 
Monte Carlo or Molecular Dynamics (MD) simulation 
free energies for both the liquid and solid phases. Re- 
cently Potekhin et al. have made accurate calculations 
of the free energy of liquid mixtures [TS],[2n] and Medin 
and Cumming calculated the phase diagram for both bi- 
nary mixtures such as C/0 and much more complicated 
multi-component mixtures |21) . 

All of these works determine liquid-solid phase equilib- 
ria by equating liquid and solid free energies that have 
been calculated separately. This procedure allows the use 
of smaller Monte Carlo or MD simulations where only a 
single phase is present at a time. However, it may be very 
sensitive to any small errors in the free energy difference 
between liquid and solid phases. Indeed for the C/0 sys- 
tem, Segretain et al. predict higher melting temperatures 
and a spindle type phase diagram while both Ogata et 
al. and Medin and Cumming predict lower melting tem- 
peratures and an azeotrope type phase diagram. In a 
spindle-type phase diagram the melting temperature of 
the mixture is always greater than the melting temper- 
ature of pure carbon, while in an azeotrope type phase 
diagram the melting temperature of the mixture can be 
lower than that of pure carbon. This difference in phase 
diagrams could be due to small errors in Segretain's solid 
free energies. 



2 



Furthermore, equating the free energies of liquid and 
solid phases provides no information on the dynamics 
of the phase transition. For example, although there 
have been some studies of nucleation for one component 
plasma systems [231 j there have been almost no stud- 
ies of nucleation for multi-component plasmas. Finally, 
interface properties are not addressed. For one compo- 
nent systems there has been some work on surface prop- 
erties, see for example For multi-component sys- 
tems there is, in general, a gradient in composition across 
the liquid-solid interface. However, the spatial extent of 
this gradient has not been determined. 

Recently, we performed direct two phase molecular dy- 
namics simulations of liquid-solid phase equilibria both 
for carbon / oxygen mixtures in WD [26j and for a com- 
plex 17 component mixture modeling the crust of an ac- 
creting neutron star [10]. These simulations have both 
liquid and solid phases present simultaneously. This al- 
lows a direct determination of the melting temperature, 
and the composition of the liquid and solid phases from a 
single simulation. Furthermore, phase equilibria for very 
complicated systems can be simulated in this way. 

However, direct simulations need to address potential 
systematic errors from finite size effects and from the 
lack of thermodynamic equilibration. Finite size effects 
are potentially important because one must fit not only 
liquid and solid phases but also two liquid-solid interfaces 
within the simulation volume. This, in general, requires 
a larger simulation volume than for simulations of only 
a single phase. However, recent computer advances have 
dramatically reduced the computational limitations on 
these larger simulations. It is now "easy" to simulate 
much larger systems than have typically been run in the 
past. 

One must run these two-phase simulations long enough 
to insure that the phases have come into thermody- 
namic equilibrium. This may require impurities to dif- 
fuse throughout the solid phase. However, diffusion in 
the solid phase is relatively fast, for these Coulomb sys- 
tems, because the ions do not have hard core interac- 
tions. There is only a relatively soft 1/r interaction be- 
tween ions. As a result, ions can move past one another. 
We have studied diffusion in Coulomb crystals in a re- 
cent paper |27| . If finite size and non-equilibrium effects 
are addressed, direct two phase simulations should yield 
accurate results. The systematic errors from two-phase 
simulations are likely very different from previous free en- 
ergy calculations. Therefore, comparing the two methods 
provides an important check on both approaches. 

In the laboratory, one can observe complex (or dusty) 
plasma crystals. Complex plasmas (CP) are low tem- 
perature plasmas containing charged micro-particles, for 
a review see Fortov et al. [28]. Often the micro- 
particles are micron sized spheres that acquire large elec- 
tric charges and the strong Coulomb interactions between 
micro-particles can lead to crystallization. Indeed plasma 
crystals were first observed in the laboratory in 1994 [25] . 
Alternatively, one can study systems of charged colloidal 



spheres. For example, Lorenz and Palberg observed melt- 
ing and freezing lines for a binary mixture of colloidal 
spheres [5U] . 

In this paper we study freezing of binary mixtures of 
both C/0 and 0/Se. The C/0 system is important for 
White Dwarfs while the 0/Se system provides a simple 
binary model of the complex rp ash composition in ac- 
creting NS. We perform MD simulations with both 27648 
and 55296 ions. This allows us to study finite size effects. 
We discuss our MD formalism in Section |ll] present re- 
sults in Section [ITT] for the carbon-oxygen, and Section [TVl 
for the oxygen-selenium system, and conclude in Section 
lYl 



II. FORMALISM 



In Section |II A| we describe our two-phase MD simu- 
lation formalism. This is very similar to what we used 
earlier for the freezing of rapid proton capture nucleosyn- 
thesis ash on accreting neutron stars (lOj and for carbon 
and oxygen mixtures in WD [25] . Next in Section II B 



we 



describe our analysis procedure for determining if a given 
region of the simulation is in a liquid or solid phase. 



A. MD formalism 

We consider a system of ions, of two different charges, 
and electrons. The electrons are assumed to form a de- 
generate Fermi gas. The ions are fully pressure ionized 
and interact with each other via screened Coulomb in- 
teractions. The potential between the ith and jth ion is 
assumed to be, 



-r/A 



(1) 



Here the ion charges are Zi and Zj , r is their separation 
and the electron screening length is A. For cold rela- 
tivistic electrons, the Thomas Fermi screening length is 
— 2a^^'^kF /t^^^^ where the electron Fermi momen- 
tum kp is kp = (STT^rie)^/^ and a is the fine structure 
constant. Finally the electron density Ue is equal to the 
ion charge density, Ug = {Z)n, where n is the ion density 
and (Z) is the average charge. We cutoff the potential 
for r > 8A. Our simulations are classical and we have 
neglected the electron mass (extreme relativistic limit). 
This is to be consistent with our previous work on neu- 
tron stars. However, the electron mass is important at 
the lower densities in WD and this may change our re- 
sults slightly [jr. Also quantum effects could play some 
role at high densities [32 ,[33|. 

The simulations can be characterized by an average 
Coulomb parameter F, 



(Z5/3)e 
a„T 



(2) 
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Here {Z^^^) is an average over the ion charges, T is the 
temperature, and the electron sphere radius Oe is = 
(3/47rn,)i/3. 

Time can be measured in units of one over the plasma 
frequency cOp . Long wavelength fluctuations in the charge 
density can undergo oscillations at the plasma frequency. 
This depends on the ion charge Z and mass M. For 
mixtures we define a hydrodynamical plasma frequency 
ijp from the simple averages of Z and M, 



Ane^iZyn 



(M) 



1/2 



(3) 



Note that other choices for the average over composition 
in Eq. [3] are possible. However, they are expected to give 
very similar results for the average plasma frequency. 

All of our carbon-oxygen simulations are run for the 
same electron density of — 5.026 x 10"'' fm~'^, while 
the oxygen-selenium simulations are run for rie ~ 2.254 x 
10^"^ fm""^. Since the pressure is dominated by the 
electronic contribution, constant electron density corre- 
sponds, approximately, to constant pressure. The density 
can be scaled to any desired value by also changing the 
temperature T so that the value of F, see Eq. [2j remains 
the same. 



B. Interface finding algorithm 

In this Section we describe an algorithm for specifying 
if a given region of the simulation belongs to the bulk 
liquid, or bulk solid, phase or if the region belongs to a 
liquid-solid interface. Often determining whether a clus- 
ter of ions is a liquid or a bcc solid is simple when visually 
inspected. However this determination is difficult to ob- 
tain numerically. For an entire system, phase determina- 
tion can be accomplished by computing the global order 
parameter Qe [34j, see below. In this work, we use the 
prescription laid out by ref. ^Sj to determine whether 
individual ions are liquid-like or solid-like. 

For each ion i, an ion j is defined as neighbors if it 
is within a given radius rcut = 4a with a the ion sphere 
radius a = (3/4?™)^/^. The vectors joining neigh- 
bors are called bonds. The direction of these vectors can 
be described by 9ij and (jjij in the frame of ion i. The 
local structure around ion i can be characterized using 
spherical harmonics Yim{dij , (jjij) by 



qimii) 
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where Nh{i) is the number of ions bonded with ion i and 
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if rij < 4a and a{rij) 



4a 



otherwise. 



(5) 



These local order parameters are large in both the solid 
and the liquid. The global order parameter Qq is calcu- 
lated from an average over all of the A'^ ions, 
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This is large in the solid due to the fact that the qQm{i) 
add up coherently. In the liquid, q^mii) add incoherently, 
so Qq is small. This coherence is exploited to determine 
local order. 

For each q^mii) a normalization is applied, 



g6m(«) 



96m (i) 



1/2 



(8) 



A dot product can now be defined of the vectors qg for 
neighboring particles i and 



(9) 



By construction, qg(i) ■ qgl*) — 1- 

We use the same criterion as ref. |35j for deter- 
mining whether two particles are connected, namely 
qg(j) • qg(j) > 0.5. This criterion will be met for most 
of the bonds in the solid. In the liquid, two neighbors 
may be in phase and considered connected, but that is 
unlikely to be true for all of the neighbors. Therefore, we 
use a threshold on the number of connections to deter- 
mine if an ion is solid-like or liquid-like. This threshold 
is 20 connections. An ion is considered solid-like if it is 
connected to 20 or more neighbors. Otherwise the ion is 
considered liquid-like. Note that on average an ion in the 
carbon-oxygen system has about 62 neighbors. 

Now that each ion is tagged as either solid-like or 
liquid-like, the interface in our two-phase simulations can 
be found. Deep in the solid, a vast majority of the ions 
within a certain radius of a given ion are identified as 
solid-like. In the bulk of the liquid a similar majority is 
identified as liquid-like. Along the interface, there is a 
mixture of solid-like and liquid-like ions. For this reason, 
we tag an ion as being in the solid or liquid if a large 
majority (> 80%) of the ions within 4a of a given ion 
are the same phase. If this criterion is not met, then the 
ion is determined to be in the interface. Performing this 
identification leads to results that will be shown in Fig. 
[2] of Section Hill Notice that ions determined to be in the 
interface are found where one would expect them, along 
the border of the solid and liquid phases. 

Note that this procedure is slightly modified for the 
oxygen-selenium system. Bonds involving oxygen ions, in 
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a predominately selenium crystal, are not as well ordered 
as selenium-selenium bonds. Therefore, for this system 
we only consider selenium-selenium bonds, as discussed 
in Sec. |lVl 



III. RESULTS FOR CARBON AND OXYGEN 
SYSTEM 

In this section we present results for the phase diagram 
and diffusion constants for carbon and oxygen systems. 
Our previous carbon and oxygen results in ref. were 
based on MD simulations with 27648 ions. Here we per- 
form three larger simulations with 55296 ions in order to 
study finite size effects. In addition we calculate diffusion 
constants in order to monitor non-equilibrium effects. 

To minimize finite size effects, we use a rectangular 
simulation volume that is twice as long in the z direc- 
tion compared to the x or y directions. We use periodic 
boundary conditions in this rectangular box. Note that, 
we evaluate the interaction between two particles as the 
single interaction with the nearest periodic image. We do 
not include an Ewald sum over further periodic images 
because our box is so large that interactions with periodic 
images other than the nearest one are very small. The 
initial conditions consist of a cube of crystalline phase 
that is stacked together in the z direction with an equal 
sized cube of liquid phase. This rectangular geometry 
increases the distance between the two liquid-solid inter- 
faces compared to a cubical simulation volume. 



A. Run with 75% oxygen 

Our first simulation has an average composition of 75% 
oxygen and 25% carbon (by number). We independently 
prepare solid and liquid initial conditions and then com- 
bine them to obtain the full 55296 ion initial conditions. 
Based on our earlier results suggesting that the solid 
phase will be enriched in oxygen from |26j . we prepare 
the liquid configuration with 70% oxygen, 30% carbon 
and the solid configuration with 80% oxygen, 20% car- 
bon. We start with a 3456 ion configuration with random 
positions and velocities and cool the system until it so- 
lidifies. We then combine 8 copies of this 3456 ion solid 
to make a 27648 ion solid configuration and continue to 
evolve this 27648 ion system for a time tojp = 8900 at 
r = 213.1. We form a liquid configuration by starting 
with an independent 3456 ion system with random coor- 
dinates and evolve the system for a time tuip — 4400 at 
the same F. We combine 8 copies of this liquid configu- 
ration to make a 27648 ion liquid and evolve for a further 
tijp — 8800. Finally we combine the 27648 ion solid con- 
figuration with the 27648 ion liquid configuration to form 
the full 55296 ion initial condition. 

This initial condition is not equilibrated for a num- 
ber of reasons. First the two interfaces between liquid 
and solid may have high energies because we have sim- 




FIG. 1: (Color on line) Final configuration of carbon ions 
(light red) and oxygen ions (black) in a 55296 ion simulation 
that consisted of 75% oxygen and 25% carbon. 



ply combined two different initial conditions. There may 
be liquid and solid ions close to each other. In addition, 
the solid part of the initial condition may not be equili- 
brated because it has carbon and oxygen ions positioned 
on lattice sites more or less at random. The equilibrated 
solid may have important correlations between ions of 
different charges. Finally, the compositions of the liquid 
and solid phases may be wrong. Therefore, carbon ions 
may diffuse into or out of the solid region until the com- 
position of the liquid and solid phases equilibrate. This 
may require evolving the system for a considerable time. 
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FIG. 2; (Color on line) Final configuration of liquid (dark 
purple), interface (light red), and solid (black in bottom half) 
ions from a 55296 ion simulation that consisted of 75% oxygen 
and 25% carbon. These liquid, solid and interface, regions 
are determined using a bond angle metric as discussed in Sec. 

mm 



We evolved the full 55296 ion system for a total time 
2.8 X 10^/a)p using a velocity Verlet "36] time step of 
0.177/cjp. We employed a simple hybrid Open MP / MPI 
computer code and use about 768 cores on the Cray XT5 
system Kraken [37] . This evolution took a total of about 
2 weeks. During the run we rescaled the velocities every 
100 time steps in order to keep the temperature approx- 
imately constant. We slowly adjust the temperature so 



that approximately half of the system remains liquid and 
half solid. In Fig. [ijwe show the final configuration of the 
55296 ions. The bottom half of the simulation volume is 
seen to contain a crystalline region. We use the proce- 
dure of Sec. II B to determine which parts of the system 
are liquid, solid, or belong to the two liquid-solid inter- 
faces. This is shown in Fig. [2j Note that the interface 
regions are not rectangular and show some fiuctuations. 

In Fig. [3] we show the fluctuations in F (or equivalently 
one over the temperature) during the run. We calculate 
the final f value as the average of F over the last third 
of the run, see Table |l] We define the scaled fluctuations 
in F as <5F, 



ST = 



F-F 



(10) 



As the run started, it was necessary to use a low tempera- 
ture (high F) in order to keep the badly non-equilibrated 
solid frozen. However as the system rapidly equilibrated 
the temperature could be quickly reduced towards its fi- 
nal equilibrated value. The temperature is then seen to 
fiuctuate around this value at the half percent level or 
less. Fig. [3] also shows the fraction of the system that 
is solid, liquid, or interface vs time. The interface frac- 
tion is very constant (since the thickness of the interface 
depends on our definition in Sec. |II B| but is time in- 
dependent). The temperature is adjusted to keep the 
fraction liquid and fraction solid more or less constant. 



XO 


f 


Xq Xq Xq 


0.75 
0.50 
0.25 


204.5(8) 
221.7(9) 
213.2(7) 


0.806(1) 0.699(3) 0.740(7) 
0.552(1) 0.454(3) 0.494(8) 
0.250(2) 0.249(2) 0.252(6) 



TABLE I: Equilibrium compositions of our 55296 ion runs. 
Oxygen number fraction of the whole system is xo, Coulomb 
parameter averaged over the last third of the run F (deter- 
mined from the final temperature and xo)- The composition 
of the solid is xq , the liquid is Xq , and the interface regions 
is Xq. Statistical errors in the last digit are quoted in paren- 
theses. 

In Fig. m we show the composition (oxygen number 
fraction xq) of the solid Xq , liquid Xq and interface Xq 
regions as a function of time. The interface composition 
is seen to remain near the average value Xq « 0.75 while 
the solid becomes enriched in oxygen so that Xq ss 0.8 
and the liquid is depleted in oxygen Xq ~ 0.7. The com- 
positions, averaged over the final third of the run, are 
collected in Table HI 

In Fig. [5] we show the oxygen composition xq vs posi- 
tion in the simulation volume at the end of the run. We 
divide up the simulation volume into fifty slices accord- 
ing to the z coordinate, with slice 1 being at the bottom 
of Fig. [1] and slice 50 being at the top. For reference 
we show in Fig. [s] the fraction of ions in a given slice 
that are in the solid, liquid, and interface regions. For 
example, we see that slices 6 to 21 are all solid. We note 
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FIG. 3: (Color on line) Fractional fluctuations in Coulomb pa- 
rameter ST (solid black line and right hand scale) , see Eq. [lO] 
vs simulation time t for a 55296 ion system with 75% oxygen 
and 25% carbon. Also shown are the fraction of the system 
that is solid (red squares), liquid (blue circles) or interface 
(brown triangles). 
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FIG. 5: (Color on line) Number fraction of oxygen xo for 
difl^erent slices (regions, see text ) of the simulation volume 
(middle panel). Also show in the top panel are diffusion con- 
stants D for carbon ions (open squares) and for oxygen ions 
(filled squares). Finally the lower panel shows the fraction 
of ions that are solid (squares), liquid (circles), or interface 
(triangles). This is averaged over the last third of the run of 
the 55296 ion simulation that is overall 75% oxygen and 25% 
carbon. 



FIG. 4: (Color on line) Number fraction of oxygen vs simula- 
tion time for liquid xq (blue circles), solid Xq (red squares) 
and interface Xq (brown triangles) regions. This is for a 55296 
ion simulation that is overall 75% oxygen and 25% carbon. 



that there are some fluctuations in the composition of 
the solid region, as a function of position, and there are 
gradients in composition across the interface regions. 

We calculated diffusion constants using the methods of 
ref. [57] in order to check on equilibration. Fig. [s] shows 
diffusion constants compared to a reference value Dq, 



for both carbon and oxygen ions at various positions in 
the simulation volume, see also Table |llj Diffusion is seen 
to be relatively fast in the liquid region vifith Z?^ for car- 
bon being only slightly larger than Dq for oxygen ions. 
This is consistent with the results of ref. [38]. In con- 
trast in the solid is almost 100 times smaller than 
in the liquid. This is similar to the one component solid 
diffusion results of ref. |27| . Furthermore in the solid, 



Dq is much smaller than D^. The composition of the 
solid can equilibrate by carbon ions diffusing in to re- 
duce Xq or diffusing out to increase Xq. One may not 
need to wait for the oxygen ions to diffuse through out 
the solid. Therefore, we expect the equilibration time of 
our system to be determined by Dq in the solid instead 
of the smaller Dq. We find that diffusion is isotropic in 
the interior of the liquid and solid regions and somewhat 
non-isotropic near the interfaces. 



Xo 






0.75 
0.50 
0.25 


0.80(1) 0.64(1) 
0.67(1) 0.55(1) 
0.64(1) 0.52(1) 


0.011(3) 0.0017(1) 
0.0070(4) 0.0012(1) 
0.0031(2) 0.0007(1) 



TABLE II: Diffusion coefficients of liquid and solid phases 
averaged over the last third of the runs. Results are expressed 
as in units of Do- The letter p denotes the phase, s for 
solid and I for liquid, while X stands for ion species, C for 
carbon and O for oxygen. Statistical errors in the last digit 
are quoted in parentheses. 



B. Run with 50% oxygen 

Our second 55296 ion run has an overall composition 
of 50% oxygen and 50% carbon. We started it in a very 
similar manner as for the 75% oxygen run except that the 
initial composition of the solid was assumed to be 55% 
oxygen, while the liquid was assumed to be 45% oxygen. 
This is based on our earlier results |26|. Figure [6] shows 
the fluctuations in F, and the fractions of the system that 
are liquid, solid, and interface as a function of time. The 
compositions of these regions vs time are shown in Fig. [7] 
and the compositions averaged over the final third of the 
run are collected in Table |lj We note that there is very 
little time dependence to these compositions. 
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FIG. 6: (Color on line) Fractional fluctuations in Coulomb pa- 
rameter ST (solid black line and right hand scale) , see Eq. 10 
vs simulation time t for a 55296 ion system with 50% oxygen 
and 50% carbon. Also shown are the fraction of the system 
that is solid (red squares), liquid (blue circles) or interface 
(brown triangles). 
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FIG. 7: (Color on line) Number fraction of oxygen vs simula- 
tion time for liquid xq (blue circles), solid Xq (red squares) 
and interface Xq (brown triangles) regions. This is for a 55296 
ion simulation that is overall 50% oxygen and 50% carbon. 



In Fig. [8]we show composition vs position. Again there 
are some fluctuations in the composition of the solid. Dif- 
fusion constants D are collected in Table [HI In general D 
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FIG. 8: Color on line) Number fraction of oxygen xo for 
different slices (regions, see text ) of the simulation volume 
(middle panel). Also show in the top panel are diffusion con- 
stants D for carbon ions (open squares) and for oxygen ions 
(filled squares). Finally the lower panel shows the fraction 
of ions that are solid (squares), liquid (circles), or interface 
(triangles). This is averaged over the last third of the run of 
the 55296 ion simulation that is overall 50% oxygen and 50% 
carbon. 



as a function of position is very similar to our results from 
the 75% oxygen simulation. However we note that D in 
the solid is somewhat smaller than for the 75% oxygen 



C. Run with 25% oxygen 

Our final 55296 ion run has an overall composition of 
25% oxygen and 75% carbon. We started it in a very 
similar manner as for the 75% oxygen run except that 
the initial composition of the solid was assumed to be 
equal to that of the liquid, with both at 25% oxygen. 
This is because our earlier simulations found only small 
chemical separation [26 . Figure [9] shows the fluctuations 
in F, and the fractions of the system that are liquid, 
solid, and interface as a function of time. The compo- 
sitions of these regions vs time are shown in Fig. [lO] 
and the compositions averaged over the final third of the 
run are collected in Table |T] We note that all composi- 
tions are near 25% oxygen. However, there is a small 
tendency for the interface regions to be slightly enriched 
in oxygen compared to both the liquid and solid regions. 
Note that the fluctuations in the interface composition 
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are larger than the fluctuations in the compositions of 
the liquid and solid regions because the interface con- 
tains fewer ions. The composition of the liquid is very 
close to that of the solid although, at some times there 
are small fluctuations in these compositions. Of course in 
the thermodynamic limit the composition of the interface 
is not relevant. Nevertheless, our MD simulation could 
be showing a real effect where the interface is enriched 
in oxygen even if the liquid and solid have nearly equal 
compositions. 
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FIG. 9: (Color on line) Fractional fluctuations in Coulomb pa- 
rameter 5T (solid black line and right hand scale), see Eq. [lO] 
vs simulation time t for a 55296 ion system with 25% oxygen 
and 75% carbon. Also shown are the fraction of the system 
that is solid (red squares), liquid (blue circles) or interface 
(brown triangles). 
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FIG. 11: Color on line) Number fraction of oxygen xq for 
different slices (regions, see text ) of the simulation volume 
(middle panel). Also show in the top panel are diffusion con- 
stants D for carbon ions (open squares) and for oxygen ions 
(filled squares). Finally the lower panel shows the fraction 
of ions that are solid (squares), liquid (circles), or interface 
(triangles). This is averaged over the last third of the run of 
the 55296 ion simulation that is overall 25% oxygen and 75% 
carbon. 
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FIG. 10: (Color on line) Number fraction of oxygen vs simu- 
lation time for liquid x^q (blue circles), solid Xq (red squares) 
and interface x\) (brown triangles) regions. This is for a 55296 
ion simulation that is overall 25% oxygen and 75% carbon. 



In Fig. [TT]we show composition vs position. The com- 
position is nearly constant independent of position. How- 
ever, again there are some fluctuations in the composition 
of the solid. Diffusion constants D are collected in Table 
|TTj Now D in the solid, for both carbon and oxygen, are 
smaller than in the 75% or 50% oxygen runs. 



D. Carbon - Oxygen Phase diagram 



We now present the liquid-solid phase diagram implied 



by the simulations in Sections III A III B and III C Fig- 
ure [12] shows the phase diagram as a function of xq- 
The y axis is the melting temperature T divided by the 
melting temperature Tc for pure carbon. We assume the 
pure carbon system melts at F™ — 178.4 [26]. This differs 
slightly from the one component plasma result because 
we include the effects of electron screening. 

The points plotted in Fig. [12] are listed in Table [T] 
Overall our 55296 ion results are close to our previous 
results that used 27648 ion simulations [M]. However, 
there are some noticeable differences. The 25% oxygen 
simulation with 55296 ions has nearly equal liquid and 
solid compositions, while the previous 27648 ion simula- 
tion found the solid slightly enriched in oxygen and found 
a slightly lower melting temperature. The 50% oxygen 
simulation with 55296 ions has a slightly larger differ- 
ence in composition between the liquid and solid than 
previous 27648 ion simulations. Finally the 75% oxygen 
simulation with 55296 ions has a slightly higher melting 
temperature than the previous 27648 ion result. This 
suggests that finite size effects, while not zero for the 
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FIG. 12: Carbon-oxygen phase diagram plotting the compo- 
sition of the liquid phase (upper red curve or circles) that 
is in equilibrium with the solid phase (lower blue curve or 
squares). Our present results from 55296 ion simulations are 
filled symbols while the open symbols are our previous results 
with 27648 ions from ref. |26j . The curves are the results of 
Medin and Gumming |21| . 



27648 ion simulations, are relatively small. 

The agreement between our 55296 ion results and the 
model of Medin and Gumming is excellent. All three 
liquid compositions, and two of the three solid compo- 
sitions, from Table [l] are very close to Medin and Cum- 
ming's results. The only small difFerencc is that our solid 
composition from the 75% oxygen simulation is not quite 
as oxygen rich as Medin and Gumming. This one differ- 
ence could be due to a small remaining systematic er- 
ror in our simulation, or to a small error in the model 
liquid and solid free energies used by Medin and Gum- 
ming. In any case, the small finite size corrections in 
going from 27648 ion to 55296 ion simulations improve 
the agreement between our MD simulations and Medin 
and Gumming. 

What are the nature of the errors from finite size and 
or non-equilibrium effects? Although small, the errors 
appear to go in different directions for our three simula- 
tions. For the 25% oxygen run the 55296 ion simulation 
has a higher melting temperature and much smaller dif- 
ference between the composition of the liquid and solid 
compared to 27648 ion results. This is in a region of the 
phase diagram where the melting temperature is almost 
independent of composition. Therefore the equilibrium 
compositions may be very sensitive to any small errors. 

For 50% oxygen, the 55296 ion simulation has a larger 
difference between the compositions of the liquid and 
solid compared to 27648 ion results. Perhaps the simplest 
finite size effect would arise if the composition gradient 
across an interface extended over a distance comparable 
to the box size. In this case simulations with small boxes 
might have more nearly equal liquid and solid composi- 



tions than simulations with larger boxes. For example 
in Fig. [9] the gradient in composition extends over a 
distance up to perhaps as many as 10 slices and the dis- 
tance between the two liquid-solid interfaces is 25 slices. 
In contrast, the 27648 ion runs from ref. are in a 
cubical box where the distance between the two liquid- 
solid interfaces is a factor of two smaller (equivalent to of 
order 12.5 slices). Thus, there could be some small finite 
size effects for 27648 ion runs. 

Finally, there could be a statistical component to the 
errors coming from a variety of non-equilibrium effects. 
For example fiuctuations in the location of an interface 
could create a new solid region and the composition of 
this region might not have time to equilibrate. Alter- 
natively there could be composition changes from large 
fiuctuations. One could test for a variety of errors of this 
type by simply repeating these simulations a number of 
times with different initial conditions. Unfortunately, we 
have not had time to do this for the present paper. How- 
ever, we plan to do this in the future. 

The good agreement between our phase diagram and 
that of Medin and Gumming strongly suggest that the 
remaining errors in our direct MD simulation approach 
are small. And in addition, it suggests that the model 
free energies employed by Medin and Gumming are good, 
at least for the carbon-oxygen system. Furthermore, we 
conclude that the phase diagram for the carbon-oxygen 
system is accurately known. We emphasize that our di- 
rect MD simulations only work at all because diffusion 
in the solid phase is relatively fast [27]. Had diffusion 
been slow then it would be very difficult to equilibrate 
the solid phases. 

We find direct two-phase MD simulations can accu- 
rately determine liquid-solid phase equilibria. This re- 
sult is very useful because direct MD simulation can be 
applied to many other systems, including very complex 
ones. Furthermore direct MD simulations for a few com- 
positions may provide very helpful benchmarks for sim- 
pler models. Note that simulations with a somewhat 
large number of particles may be necessary and these 
may have to be run for extended simulation times. How- 
ever, rapid advances in computer power should make such 
simulations even easier in the future. 



IV. RESULTS FOR OXYGEN AND SELENIUM 
SYSTEMS 

In this section we present results for the phase dia- 
gram and diffusion constants for oxygen and selenium 
systems. This two-component system has a much larger 
ratio of charges, and this leads to a very different phase 
diagram, than for the carbon and oxygen system of Sec. 
[ml We perform MD simulations with both 27648 and 
55296 ions in order to monitor finite size effects. In gen- 
eral we present figures for the larger 55296 ion simula- 
tions and then include both 27648 and 55296 ion results 
in tables. These simulations follow closely the formalism 
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of Sec. [TTjand the procedures of Sec. Ill 

However the interface finding algorithm in Sec. |IIB| 
is slightly modified. The bond-angles for oxygen ions 
can be liquid-like (more random) even in the solid phase. 
Therefore, we only take into account selenium ions in the 
bond angle algorithm. By ignoring all oxygen ions in the 
system we determine whether a selenium ion is solid-like 
or liquid-like by analyzing how it bonds to its neighboring 
selenium ions. As in Sec. |IIB ions are considered neigh- 



bors if they are within a distance r^j < 4a of each other. 
Once all selenium ions have been identified as solid-like or 
liquid-like we tag an ion (oxygen or selenium) as being in 
the bulk of the solid or liquid if a large majority (> 80%) 
of the selenium ions within 4a of it are the same phase. 
If this criterion is not met, then the ion is determined to 
be in the interface. 

We now discuss runs with 98%, 90%, 80%, 70%, 60%, 
and 50% selenium and then we will collect results for 
the oxygen-selenium phase diagram. Although most runs 
appear to be equilibrated by the end of the simulations, 
we note that the 60% and 50% selenium runs are not 
equilibrated after using a reasonable amount of computer 
time. We discuss this below. 
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FIG. 13: (Color on line) Volume fractions of solid (squares), 
liquid (circles), and interface (triangles) versus time for a 
55296 ion simulation that is overall 2% oxygen and 98% se- 
lenium. Also shown are fluctuations SF in the Coulomb pa- 



rameter (solid line), see Eq. 10 



A. Run with 98% Selenium 

As for the carbon-oxygen systems we prepare solid 
and liquid initial conditions separately and then combine 
them to obtain the full 55296 ion initial conditions. 

To prepare the solid we start with a 432 ion system that 
is composed of 98% selenium ions with random positions 
and velocities. The system is evolved at a temperature 
close to its expected melting temperature for tuip = 4500 
using a time step of O.llwp. Due to finite size effects it 
solidifies. We then make 8 copies of this 432 ion system 
to obtain a 3456 ion solid and evolved it for tuij, ~ 20000 
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FIG. 14: (Color on line) Number fraction of selenium xse in 
the solid (squares), liquid (circles), and interface (triangles) 
versus time for a 55296 ion simulation that is overall 2% oxy- 
gen and 98% selenium. 



keeping a time step of O.lluip. Eight copies of this system 
are made to obtain a 27648 ion solid. This solid is evolved 
for tuip ~ 20000 with a time step of 0.22wp. 

The liquid is prepared in a similar fashion. We start 
with a 3456 ion system that is 98% selenium ions with 
random positions and velocities. This system is evolved 
for tuip ~ 20000 with a time step of O.lluip. Eight copies 
of this system are made to obtain a 27648 ion liquid which 
is evolved for tuip ~ 20000 using a time step of 0.22(1;^. 
Finally, we place the 27648 ion liquid configuration on 
top of the 27648 ion solid configuration to form the full 
55296 ion initial conditions. 

The temperature was adjusted during the run to keep 
approximately equal volume fractions of solid and liquid 
as shown in Fig. 
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The number fraction of selenium 
xse in the solid, liquid and interface are shown in Fig. 
[M] versus simulation time, see also Table [Till These frac- 
tions do not change much in the second half of the run 
suggesting that the system is (at least approximately) 
equilibrated. 





r 






4e 


0.98 


201.1(6) 


0.9926(5) 


0.966(1) 


0.975(2) 


0.90 


213.4(8) 


0.9708(7) 


0.843(1) 


0.884(4) 


0.80 


252.5(6) 


0.968(1) 


0.681(1) 


0.762(7) 


0.70 


289(1) 


0.970(1) 


0.547(1) 


0.655(4) 


0.60 


391(1) 


0.935(2) 


0.408(1) 


0.513(8) 


0.50 


459(2) 


0.902(2) 


0.308(1) 


0.436(9) 



TABLE III: Equilibrium compositions of our 55296 ion 
runs. Selenium number fraction of the whole system is xse, 
Coulomb parameter averaged over the last third of the run F 
(determined from the final temperature and xo). The com- 
position of the solid is Xs, the liquid is xi, and the interface 
regions is Xi. Statistical errors are quoted in the last digit in 
parentheses. 



In Fig. 15 in the middle panel, we show the number 
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fraction of selenium xse for different slices (regions) of 
the simulation volume. This is at the end of the sim- 
ulation. In the solid, xse is seen to be nearly constant 
and independent of position. This is consistent with the 
system having reached equilibration. Diffusion constants 
D are also shown in Fig. [15] in the top panel and listed in 
Table ITVl For selenium D in the solid is seen to be three 
orders of magnitude smaller than D in the liquid. How- 
ever the behavior is very different for oxygen. Diffusion 
is nearly the same in the liquid and solid regions! Indeed 
in the solid, D for oxygen is over three orders of magni- 
tude larger than D for selenium. Presumably the effec- 
tive size of oxygen ions in the solid (ion sphere radius) 
is small enough so that the oxygen can diffuse relatively 
easily through the larger crystal lattice of selenium ions. 
Note that this behavior for oxygen in a selenium crystal 
is very different from that for carbon in an oxygen crystal 
as found in Sec. |IIH see Fig. [5] 
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FIG. 15: (Color on line) Diffusion constant D (top panel)) 
for selenium (open squares) and oxygen ions (filled squares) 
for different slices (regions) of the simulation volume. Num- 
ber fraction of selenium xse (middle panel) for different slices. 
Fraction of ions x (lower panel) that are in the solid (squares), 
liquid (circles) and interface (triangles) for different slices. 
This is an average over the last 1.4 x lO^uipt of a 55296 ion 
simulation that is overall 2% oxygen and 98% selenium. 







Dh Die 


0.98 


3.17(1) 0.736(3) 


5.66(1) 0.00080(2) 


0.90 


3.04(1) 0.729(2) 


5.34(3) 0.00064(7) 


0.80 


2.88(1) 0.721(6) 


4.50(9) 0.00044(7) 


0.70 


2.81(1) 0.771(2) 


1.82(1) 0.00014(1) 


0.60 


2.36(1) 0.617(6) 


0.18(2) 0.00002(1) 


0.50 


2.12(1) 0.547(2) 


0.1(1) 0.00001(1) 



TABLE IV: Diffusion coefficients of liquid and solid phases 
averaged over the last third of the runs. Results are expressed 
as in units of Do. The letter p denotes the phase, s for 
solid and I for liquid, while X stands for ion species, O for 
oxygen and 5*6 for Selenium. Statistical errors are quoted in 
the last digit in parentheses. 



solid and 3456 ion liquid systems were set to 90% sele- 
nium. We then follow the procedure laid out in Sec. |IV A| 
to obtain a 55296 ion configuration. 

The temperature was adjusted during the run to keep 
approximately equal volume fractions of solid and liquid 
as shown in Fig. |16| The number fraction of senium xse 
in the solid, liquid and interface are shown in Fig. [17] 
versus simulation time. These fractions do not change 
much in the second half of the run suggesting that the 
system is (at least approximately) equilibrated, as for the 
run with 98% selenium. 
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FIG. 16: (Color on line) Volume fractions of solid (squares), 
liquid (circles), and interface (triangles) versus time for a 
55296 ion simulation that is overall 10% oxygen and 90% 
selenium. Also shown are fluctuations 8V in the Coulomb 
parameter (solid line), see Eq. |10[ 



B. Run with 90% Selenium 

The inital conditions for the 90% selenium simulations 
were prepared similarly to the 98% ones. The only dif- 
ference was that the initial compositions of the 432 ion 



In Fig. 18 



in the middle panel, we show the number 
fraction of selenium xse for different slices (regions) of 
the simulation volume. This is at the end of the sim- 
ulation. In the solid, xse is seen to be nearly constant 
and independent of position. This is consistent with the 
system having reached equilibration. Diffusion constants 
D are also shown in Fig. [18] in the top panel and are seen 



to behave in a similar way to Fig. 15 
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FIG. 17: (Color on line) Number fraction of selenium xse in 
the solid (squares), liquid (circles), and interface (triangles) 
versus time for a 55296 ion simulation that is overall 10% 
oxygen and 90% selenium. 



riched in selenium and desiring the system to reach equi- 
librium faster, we start the liquid and the solid subsys- 
tems with different initial compositions. We start with 
a 432 ion solid configuration that is 90% selenium and 
a 3456 liquid configuration that is 70% selenium. We 
then follow the procedure of Sec. |IV A| to obtain the 
55296 80% selenium system. We note here that due to 
their different initial compositions, the solid and the liq- 
uid initally have different electron densities. However the 
electron densities quickly equilibrate. 

The temperature was adjusted during the run to keep 
approximately equal volume fractions of solid and liquid 
as shown in Fig. 19 The number fraction of selenium 



xse in the solid, liquid and interface are shown in Fig. 
[20]versus simulation time. These fractions do not change 
much in the second half of the run suggesting that the 
system is (at least approximately) equilibrated. 
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FIG. 18: (Color on line) Diffusion constant D (top panel)) 
for selenium (open squares) and oxygen ions (filled squares) 
for different slices (regions) of the simulation volume. Num- 
ber fraction of selenium xse (middle panel) for different slices. 
Fraction of ions x (lower panel) that are in the solid (squares), 
liquid (circles) and interface (triangles) for different slices. 
This is an average over the last 1.4 x lO^cjpt of a 55296 ion 
simulation that is overall 10% oxygen and 90% selenium. 
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FIG. 19: (Color on line) Volume fractions of solid (squares), 
liquid (circles), and interface (triangles) versus time for a 
55296 ion simulation that is overall 20% oxygen and 80% 
selenium. Also shown are fiuctuations 5V in the Coulomb 



parameter (solid line), see Eq. 10 
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C. Run with 80% Selenium 

The inital conditions for this simulation were prepared 
similarly to the initial conditions for the 98% and 90% 
selenium systems. However, expecting the solid to be en- 



FIG. 20: (Color on line) Number fraction of selenium xse in 
the solid (squares), liquid (circles), and interface (triangles) 
versus time for a 55296 ion simulation that is overall 20% 
oxygen and 80% selenium. 
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In Fig. 21 in the middle panel, we show the number 
fraction of selenium xse for different slices (regions) of 
the simulation volume. This is at the end of the sim- 
ulation. In the solid, xse is seen to be nearly constant 
and independent of position. This is consistent with the 
system having reached equilibration. Diffusion constants 
D are also shown in Fig. [21] in the top panel and are seen 
to behave in a similar way to Fig. |15| 



10" r 



D 
Da 



10-1 
10-= 
10-3 
10-* 
1.00 

0.90 [f 
0.80 \f 

0.70 «r 















r Q 


□ 






□ 






° Dse/Do 


B 


!" '^'-ttnEEci 


^o/Do 


u 




4^' 



1.00 
0.80 
0.60 
0.40 
0.20 
0.00 [Ptele 



. gpOOOOOCTOOaCPGOOQOOCTQ . 



10 



20 30 
slices 



40 



50 



FIG. 21: (Color on line) Diffusion constant D (top panel)) 
for selenium (open squares) and oxygen ions (filled squares) 
for different slices (regions) of the simulation volume. Num- 
ber fraction of selenium xse (middle panel) for different slices. 
Fraction of ions x (lower panel) that are in the solid (squares), 
liquid (circles) and interface (triangles) for different slices. 
This is an average over the last 1.4 x lO^ujpt of a 55296 ion 
simulation that is overall 20% oxygen and 80% selenium. 



However, it is possible that xse could continue to change 
for much longer run times. 
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FIG. 22: (Color on line) Volume fractions of solid (squares), 
liquid (circles), and interface (triangles) versus time for a 
55296 ion simulation that is overall 30% oxygen and 70% 
selenium. Also shown are fluctuations SF in the Coulomb 



parameter (solid line), see Eq. 10 
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FIG. 23: (Color on line) Number fraction of selenium xse in 
the solid (squares), liquid (circles), and interface (triangles) 
versus time for a 55296 ion simulation that is overall 30% 
oxygen and 70% selenium. 



D. Run with 70% Selenium 

The initial conditions for this simulation were prepared 
similarly to the initial conditions for the 80% selenium 
system. The initial 432 ion solid was, again, composed 
of 90% selenium while the 3456 ion liquid was set to 50% 
selenium. 

The temperature was adjusted during the run to keep 
approximately equal volume fractions of solid and liquid 
as shown in Fig. |22[ The number fraction of selenium 
xse in the solid, liquid and interface are shown in Fig. 
1231 versus simulation time. The fraction of selenium in 
the solid may be increasing very slowly with time while 
xse in the liquid may be slowly decreasing. Because this 
change is very slow the system may be near equilibration. 



In Fig. 24 in the middle panel, we show the number 
fraction of selenium xse for different slices (regions) of 
the simulation volume. This is at the end of the sim- 
ulation. In the solid, xse is seen to be nearly constant 
and independent of position. This is consistent with the 
system having reached equilibration. Diffusion constants 
D are also shown in Fig. [24] in the top panel and are seen 
to behave in a similar way to Fig. 



15 However now D 



for selenium is seen to change slightly with position in 
the solid near the interface regions. 



E. Run with 60% Selenium 

The initial conditions for this simulation were prepared 
similarly to the initial conditions for the 80% selenium 
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FIG. 24: (Color on line) Diffusion constant D (top panel)) 
for selenium (open squares) and oxygen ions (filled squares) 
for different slices (regions) of the simulation volume. Num- 
ber fraction of selenium xse (middle panel) for different slices. 
Fraction of ions x (lower panel) that are in the solid (squares), 
liquid (circles) and interface (triangles) for different slices. 
This is an average over the last 1.4 x IQ^uipt of a 55296 ion 
simulation that is overall 30% oxygen and 70% selenium. 
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FIG. 25: (Color on line) Volume fractions of solid (squares), 
liquid (circles), and interface (triangles) versus time for a 
55296 ion simulation that is overall 40% oxygen and 60% 
selenium. Also shown are fluctuations 5V in the Coulomb 



parameter (solid line), see Eq. 10 
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system. The initial 432 ion solid was composed of 90% 
selenium while the 3456 ion liquid was set to 30% sele- 



FIG. 26: (Color on line) Number fraction of selenium xse. in 
the solid (squares), liquid (circles), and interface (triangles) 
versus time for a 55296 ion simulation that is overall 40% 
oxygen and 60% selenium. 



The system started with a large solid fraction. The 
temperature was increased to bring the solid fraction ap- 
proximately equal with the liquid fraction as shown in 
Fig. |25[ The number fraction of senium xse in the solid, 
liquid and interface are shown in Fig. [26] versus simu- 
lation time. Again the fraction of Selenium in the solid 
may be increasing slightly over the second half of the run. 

In Fig. |27[ in the middle panel, we show the number 
fraction of selenium xg^^ for different slices (regions) of 
the simulation volume. This is at the end of the simula- 
tion. In the solid, xse depends on position, as do diffusion 
constants D as shown in the top panel of Fig. [27] This 
clearly indicates that the system has not reached equi- 
librium. Therefore we can not use this run to determine 



the phase diagram in Sec. IV G Unfortunately, the equi- 
libration time may be very long for this composition and 
therefore require an unreasonable amount of simulation 
time in order to bring the system in equilibration. 



F. Run with 50% Selenium 

The initial conditions for this simulation were prepared 
similarly to the initial conditions for the 80% selenium 
system. The initial 432 ion solid was, again, composed 
of 90% selenium while the 3456 ion liquid was set to 10% 
selenium. 

The system started with a large solid fraction. The 
temperature was increased to bring the solid fraction ap- 
proximately equal with the liquid fraction as shown in 
Fig. |28| The number fraction of selenium xse in the 
solid, liquid and interface are shown in Fig. [29] versus 
simulation time. Now the fraction of Selenium in the 
solid may be decreasing slightly, and xse for the liquid 
increasing slightly, with time over the second half of the 
run. 

In Fig. |30[ in the middle panel, we show the number 
fraction of selenium xse for different slices (regions) of 
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FIG. 27: (Color on line) Diffusion constant D (top panel)) 
for selenium (open squares) and oxygen ions (filled squares) 
for different slices (regions) of the simulation volume. Num- 
ber fraction of selenium xse (middle panel) for different slices. 
Fraction of ions x (lower panel) that are in the solid (squares), 
liquid (circles) and interface (triangles) for different slices. 
This is an average over the last 1.4 x 10^ (Dpt of a 55296 ion 
simulation that is overall 40% oxygen and 60% selenium. 



the simulation volume. This is at the end of the sim- 
ulation. In the solid, xse depends on position, as do 
diffusion constants D as shown in the top panel of Fig. 
|30| This is similar to the run with 60% selenium and 
clearly indicates that the system has not reached equilib- 
rium. Therefore we can not use this run to determine the 
phase diagram in Sec. |IV G[ Unfortunately, the equilibra- 
tion time for this composition may also be very long and 
therefore require an unreasonable amount of simulation 
time in order to bring the system in equilibration. 



G. Oxygen-Selenium Phase diagram 

We now present the liquid-solid phase diagram implied 
by th e si mulati ons in Sections [WTl [IVB| [iVCl [TVD 



rVEl and [TVFl We use the data in Table |III| Figure |31 
shows the phase diagram as a function of xse- The y 
axis is the melting temperature T divided by the melting 
temperature Tq for pure oxygen. We assume the pure 
oxygen system melts at Tm = 178.4 [26]. This differs 
slightly from the one component plasma result because 
we include the effects of electron screening. 

The filled upward pointing triangles in Fig. [31] show 
the composition of the liquid phase, and the filled squares 
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FIG. 28: (Color on line) Volume fractions of solid (squares), 
liquid (circles), and interface (triangles) versus time for a 
55296 ion simulation that is overall 50% oxygen and 50% 
selenium. Also shown are fluctuations ST in the Coulomb 



parameter (solid line), see Eq. 10 
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FIG. 29: (Color on line) Number fraction of selenium xse in 
the solid (squares), liquid (circles), and interface (triangles) 
versus time for a 55296 ion simulation that is overall 50% 
oxygen and 50% selenium. 



the composition of the solid phase for 55296 ion simula- 
tions. Also shown as filled downward pointing triangles, 
and filled diamonds are the liquid and solid compositions 
for runs that are clearly not equilibrated. These points 
should not be used in the determination of the phase di- 
agram. The open triangles (liquid) and squares (solid) 
show results for smaller 27648 ion simulations, see Table 
[V] In general there is very good agreement between equi- 
librated 55296 and 27648 ion simulations. This suggests 
that finite size effects, while not strictly zero, are small. 

The blue solid lines in Fig. [3l] show the solid com- 
position, and the dotted brown line the liquid composi- 
tion, for the phase diagram of Medin and Gumming |21) . 
There is qualitative agreement between these curves and 
our results for the xse > 0.5 half of the phase diagram 
where we have apparently equilibrated independent re- 
sults. However there are some differences in detail. We 
find a somewhat larger selenium solid composition while 
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FIG. 30: (Color on line) Diffusion constant D (top panel)) 
for selenium (open squares) and oxygen ions (filled squares) 
for different slices (regions) of the simulation volume. Num- 
ber fraction of selenium xse (middle panel) for different slices. 
Fraction of ions x (lower panel) that are in the solid (squares), 
liquid (circles) and interface (triangles) for different slices. 
This is an average over the last 1.4 x 10^ u>pt of a 55296 ion 
simulation that is overall 50% oxygen and 50% selenium. 



XSe 


f 


^Se ^'se ^Se 


0.98 


198.1(7) 


0.991(1) 0.968(1) 0.976(2) 


0.90 


212.0(8) 


0.969(2) 0.846(4) 0.888(4) 


0.80 


244(2) 


0.966(2) 0.683(3) 0.764(6) 


0.70 


279(1) 


0.970(2) 0.560(3) 0.650(8) 


0.65 


329(1) 


0.973(2) 0.469(3) 0.559(6) 


0.60 


383(1) 


0.971(2) 0.415(3) 0.487(7) 


0.50 


456(1) 


0.959(1) 0.330(1) 0.400(6) 



TABLE V: Equilibrium compositions of 27648 ion runs. Se- 
lenium number fraction of the whole system is xse- Note 
that runs with xse = 0.60 and 0.50 are not equilibrated. The 
Coulomb parameter averaged over the last third of the run is 
r. The composition of the solid is Xg^, the liquid is x^^^, and 
the interface regions is x\;^. Statistical errors are quoted in 
the last digit in parentheses. 



our melting temperature for the liquid is somewhat lower 
than Medin and Gumming. This difference in melting 
temperature may be due to electron screening effects that 
are included in our simulations and neglected in Medin 
and Cumming's free energies. 

Screening depends on the ratio of ion sphere radius a 
to screening length A, see Eq. [l] For a very relativistic 
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FIG. 31: Selenium-oxygen phase diagram plotting the com- 
position of the liquid phase (upper red curve or triangles) 
that is in equilibrium with the solid phase (lower blue curve 
or squares). Results from 55296 ion simulations are filled 
symbols while the open symbols are results with 27648 ions. 
Besides results discussed in the text there is an extra 27648 
simulation with xse = 0.65. The curves are the results of 
Medin and Gumming [21| . 



electron gas this ratio depends only on average ion charge 
(Z) and is independent of density, 



a 
A 



a5(Z): 



(12) 



For a one component Yukawa system the value of F at 
the melting point, F„i, is [31] 

TmiiZ)) w 171.8 + 42.46k2 + 3.841k*, (13) 

for K < 1.4. For pure oxygen we have 

F,„(8) = 178, (14) 

while for pure selenium, 

F„(34) = 188. (15) 

Thus pure selenium melts at a 6% higher F value than 
pure oxygen because of enhanced electron screening. In 
order to study the differences between our MD simulation 



results for T/Tq in Fig. 31 and Medin and Cumming's 



results, we rescale Medin and Cumming's melting tem- 
peratures according to 



T 



r„(8) T 

TU{Z)) To' 



(16) 



and plot the rescaled results in Fig. 32 This procedure 
ensures that the rescaled Medin and Cumming results 
will reproduce Eq. [14] for pure oxygen and Eq. [15] for 
pure selenium. In between, for a mixture of oxygen and 
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composition considered in ref. [TU] was predominately 
selenium, with only small concentrations of oxygen and 
a number of other impurities. We modeled this 17 com- 
ponent composition with the binary system of oxygen 
impurities mixed with the dominant element selenium. 
Direct MD simulations of the full rp ash composition in 
ref. [TU] Table I found the concentration of oxygen in the 
liquid phase to be six times larger than the oxygen con- 
centration of the solid phase. While we find, in the first 
two rows of Table that the oxygen concentration of 
the liquid phase, for our simplified binary mixture sim- 
ulations, to be five times larger than that in the solid 
phase. We conclude that this binary mixture model pro- 
vides a reasonable description of the freezing behavior of 
the rp ash. 



FIG. 32: Enlargement of the large xse part of the Selenium- 
oxygen phase diagram from Fig. |31[ The heavy dashed and 
solid lines are the results of Medin and Gumming [21] i while 
the light dashed lines are our rescaling of Medin and Cum- 
ming's results to approximately include the effects of electron 
screening as discussed in the text. 



selenium, we somewhat arbitrarily assume that the elec- 
tron screening correction depends only on (Z) . 

There is good agreement between our MD results in 
Fig. ^ 



32 and these rescaled Medin and Gumming results 



for simulations with 90% or 98% selenium. However for 
simulations with smaller xse: our MD results give some- 
what lower melting temperatures and have less oxygen 
in the solid phase. These differences could be due to 
finite size or non-equilibrium effects in our MD simula- 
tions. However, we only find very small differences be- 
tween 27648 and 55296 ion simulations. Furthermore, 
except for the low xse = 50% or 60% runs, we find very 
little time dependence in our simulations. Alternatively 
the differences could be due to electron screening effects 
for mixtures that are not well described by the rescaling 
in Eq. 16 Finally, the differences could be due to limita- 



tions in the liquid and solid free energies used by Medin 
and Gumming. 

Unfortunately our simulations with 50% and 60% se- 
lenium did not reach equilibrium. Therefore we are not 
able to effectively study the phase diagram for low sele- 
nium concentrations. Medin and Gumming predict re- 
gions of the phase diagram with equilibria between dif- 
ferent solid phases. Perhaps by carefully preparing ini- 
tial conditions that include two solid phases of differ- 
ent compositions, one may be able to study solid-solid 
phase equilibria with our direct MD simulation proce- 
dure. However, the small diffusion constants for selenium 
in the solid have made it difficult for us to equilibrate the 
simulations with small selenium concentrations presented 
here. Note that this region of the phase diagram, with 
small xse, may not be important in applications for neu- 
tron stars. 

The complex rapid proton capture nucleosynthesis ash 



V. SUMMARY AND CONCLUSIONS 

We have determined the liquid-solid phase diagram 
for carbon-oxygen and oxygen-selenium plasma mixtures 
using two-phase MD simulations. We identified liquid, 
solid, and interface regions in our simulations using a 
bond angle metric described in Sec. |IIB| To study fi- 
nite size effects, we performed both 27648 and 55296 ion 
simulations. To help monitor non-equilibrium effects, we 
calculated diffusion constants Di. For the carbon-oxygen 
system, we find that Dq for oxygen ions in the solid is 
much smaller than Df-, for carbon ions and that both dif- 
fusion constants are 80 or more times smaller than diffu- 
sion constants in the liquid phase. There is excellent 
agreement between our carbon-oxygen phase diagram 
and that predicted by Medin and Gumming !2lj. This 
suggests that errors from finite size and non-equilibrium 
effects are small, and that the carbon-oxygen phase dia- 
gram is now accurately known. 

The oxygen-selenium system, with a larger ratio of 
charges than carbon-oxygen, can serve as a simple two 
component model of the complex rapid proton capture 
ash composition on an accreting neutron star. We find 
that diffusion of oxygen in a predominately selenium 
crystal is remarkably fast and is comparable to diffusion 
in the liquid phase. Our MD simulations have a some- 
what lower melting temperature for the oxygen-selenium 
system than that predicted by Medin and Gumming. 
This is in part due to electron screening effects, that 
are included in our simulations and may be neglected 
by Medin and Gumming. In the future, we will present 
MD simulations of the phase diagram for the three com- 
ponent carbon-oxygen-neon system to include the effects 
of neon impurities in carbon-oxygen white dwarfs. 
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87ER40365 and by the National Science Foundation 
through TeraGrid resources provided by the National 
Institute for Gomputational Sciences under grant TG- 
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